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1  Introduction 


To  illustrate  the  how  one  may  implement  a  lattice-gas  with  long-range  in¬ 
teractions,  let  us  consider  for  simplicity  a  two-dimensional  example  with  a 
system  having  only  one  interaction  range  and  consider  only  an  attractive 
interaction.  The  more  general  case  of  multiple  interaction  ranges  with  both 
repulsive  and  attractive  interactions  and  in  higher  dimensions  will  follow 
directly.  Specifically,  w'e  discuss  our  new  molecular  dynamics  lattice-gas  al¬ 
gorithm  that  uses  eight  interaction  ranges  and  both  repulsive  and  attractive 
interactions  to  approximate  a  Lenard- Jones  intermolecular  potential. 

A  long-range  lattice-gas  has  been  implemented  on  the  MIT  cellular  au¬ 
tomata  machine  prototype,  the  CAM-8.  Consequently,  first  the  CAM-8  ar¬ 
chitecture  is  briefly  described.  Next  a  brief  description  of  what  a  lattice-gas 
automaton  is  and  explain  why  it  is  an  exactly  computable  representation  of  a 
dynamical  system  is  given.  One  of  the  principal  requirements  for  a  lattice-gas 
with  microscopic  finite-point  group  symmetry  to  give  rise  to  macroscopic  con¬ 
tinuous  rotational  symmetry  is  that  the  underlying  lattice  must  be  isotropic. 
Therefore  1  describe  what  it  means  for  a  lattice  to  be  isotropic.  Working  in 
two  dimensions  is  much  easier  than  working  in  three,  both  for  implement¬ 
ing  computer  models  and  for  describing  them.  For  this  reason  I  present  the 
long-range  lattice-gas  algorithm  in  two  dimensions  on  the  triangular  lattice. 

When  introduced  to  the  triangular  lattice-gas  model  for  the  first  time, 
one  inevitably  asks  the  following  question:  Why  does  discrete  dynamics  fail 
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to  reproduce  the  correct  continuum  hydrodynamic  limit  when  implemented 
on  a  square  lattice?  One  finds  that  four  momentum  states  are  insufficient 
by  noting  that  the  derivation  of  the  Navier-Stokes  equation  relies  on  the 
expansion  of  the  momentum  flux  density  tensor  in  terms  of  the  isotropic 
tensor  In  turn  the  tensor  could  be  expanded  in  products  of  two 
dimensional  Kronecker  deltas,  given  below  in  Eq.  (10).  For  the  square  lattice 
case,  the  lattice  vectors  are  orthogonal  and  E‘^  cannot  be  decomposed  into 
two-dimensional  Kronecker  deltas.  Instead 

■^tyfc/|.S=4  =  ^^ijkl 

where  5ijki  is  a  four  dimensional  Kronecker  delta,  illustrating  the  lack  of 
isotropy  of  the  momentum  flux  density  on  a  square  lattice-gas.  Since  five 
nearest  neighbors  are  not  space  filling,  the  next  possible  choice  is  six  or  the 
triangular  lattice.  The  simplest  discrete  dynamics  in  two  dimensions  is  known 
as  a  hexagonal  lattice-gas  or  an  FHP  lattice-gas,  after  its  originators  Uriel 
Frisch,  Brosl  Hasslacher,  and  Yves  Pomeau  [1].  Since  the  long-range  lattice- 
gas  still  retains  local  collisions,  the  simple  FHP  model  is  presented  here  for 
completeness.  Next  we  examine  the  long-range  2-body  interaction,  restricting 
ourselves  to  central-body  attractive  interactions,  for  the  sake  of  simplicity. 
Two  different  bound  states,  the  bounce-back  orbit  and  the  clockwise  orbit 
are  discussed. 

When  implementing  lattice-gas  algorithms  it  is  often  useful  trj'’  to  find 
economical  ways  of  expressing  the  collisions  or  interactions,  to  reduce  the 
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size  of  a  look-up  table  or  reduce  the  depth  of  the  logical  representation  of  the 
algorithm.  To  this  end  I  briefly  discuss  some  symmetries  inherent  in  long- 
range  interactions,  in  particular  I  introduce  parity  and  conjugation  symme¬ 
tries.  Finally,  I  discuss  the  implementation  of  a  multi-long-range  lattice-gas. 
Remarkably,  this  methodology  allows  us  to  model  solid-state  dynamics  and 
as  such  offers  an  alternative  to  the  usual  method  of  computational  molecular 
dynamics. 

2  The  Cellular  Automata  Machine  CAM-8 

The  cellular  automata  machine  CAM-8  architecture  devised  by  Norman 
Margolus  of  the  MIT  Laboratory  for  Computer  Science  [2,  3]  is  the  latest  in  a 
line  of  cellular  automata  machines  developed  by  the  Information  Mechanics 
Group  at  MIT  [4,  5,  6].  It  is  optimized  for  performing  lattice-gas  simulations. 
The  CAM-8  architecture  itself  is  a  simple  abstraction  of  lattice  gas  dynamics. 
Lattice  gas  data  streaming  and  collisions  are  directly  implemented  in  the 
architecture.  The  communication  network  is  a  cartesian  three  dimensional 
mesh.  Crystallographic  lattice  geometries  can  be  directly  embedded  into  the 
CAM-8.  Each  site  of  the  lattice  has  a  certain  number  of  bits  (a  multiple  of 
16)  which  we  refer  to  as  a  “cell”.  Each  bit  of  the  cell,  or  equivalently  each 
bit  plane  of  the  lattice,  can  be  translated  through  the  lattice  in  any  arbitrary 
direction.  The  translation  vectors  for  the  bit  planes  are  termed  “kicks” .  The 
specification  of  the  x,y,  and  z  components  of  the  kicks  for  each  bit  plane  (or 
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Figure  1:  MIT  Laboratory  for  Computer  Science  Cellular  Automata  Machine 
CAM-8.  This  8  module  prototype  can  evolve  a  D-dimensional  cellular  space 
with  32  million  sites  where  each  site  has  16  bits  of  data  with  a  site  update 
rate  of  200  million  per  second. 

hyperplane)  exactly  defines  the  lattice.  The  kicks  can  be  changed  during  the 
simulation.  Thus,  the  data  movement  in  the  CAM-8  is  general.  Once  the 
kicks  are  specified,  the  coding  of  the  lattice-gas  streaming  is  completed.  In 
effect,  the  kicks  determine  all  the  global  permutations  of  the  data. 

Every  configuration  of  the  local  data  within  a  cell  is  an  element  of  an 
equivalence  class  having  a  particular  value  of  the  conserved  quantities  of 
the  dynamics.  Local  permutations  of  data  in  an  equivalence  class  occur 
w'ithin  the  cells.  These  local  permutations  are  the  computational  metaphor 
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for  physical  conservative  collisions  between  particles.^  All  local  permuta¬ 
tions  are  implemented  in  look-up  tables.  That  is,  all  possible  physical  events 
with  a  certain  input  configuration  and  a  certain  output  configuration  are 
precomputed  and  stored  in  SRAM,  for  fast  table  look-up.  The  width  of  the 
CAM-8  look-up  tables  are  limited  to  16-bits,  or  64K  entries.  This  is  a  rea¬ 
sonable  width,  satisfying  the  opposing  considerations  of  model  complexity 
versus  memory  size  limitations  for  the  SRAM.  Site  permutations  of  data 
wider  than  16-bits  must  be  implemented  in  several  successiv'-e  table  look-up 
passes.  Since  the  look-up  tables  are  double  buffered,  a  scan  of  the  space  can 
be  performed  while  a  new  look-up  table  is  loaded  for  the  next  scan. 

Figure  2  is  a  schematic  diagram  of  a  CAM-8  system.  On  the  left  is  a 
single  hardware  module — the  elementary  ‘"chunk”  of  the  architecture.  On  the 
right  is  an  indefinitely  extendable  array  of  modules  (drawn  for  convenience 
as  two-dimensional,  the  array  is  normally  three-dimensional).  A  uniform 
spatial  calculation  is  divided  up  evenly  among  these  modules,  with  each 
module  able  to  simulate  a  volume  containing  millions  of  fine-grained  spatial 
sites  in  a  sequential  fashion.  In  the  diagram,  the  solid  lines  between  modules 
indicate  a  local  mesh  interconnection.  These  wires  are  used  for  spatial  data 
movements.  There  is  also  a  tree  network  (not  shown)  connecting  all  modules 
to  the  front-end  host,  a  SPARCstation  with  a  custom  SBus  interlace  card. 

^The  CAM-8  is  not  limited  to  performing  only  local  permutations;  it  can  do  general 
mappings.  However,  since  we  are  interested  only  in  particle  conserving  reversible  dynaimcs, 
local  permutations  are  sufficient. 


Figure  2:  CAM-8  System  Diagram,  (a)  A  single  processing  node,  with 
DRAM  site  data  flowing  through  a  SRAM  lookup  table  and  back  into  DRAM, 
(b)  Spatial  array  of  CAM-8  nodes,  with  nearest-neighbor  (mesh)  interconnect 
(1  wire/bit-slice  in  each  direction). 

This  SPARCstation  controls  the  CAM-8.  It  downloads  a  bit-mapped  pattern 
as  the  initial  condition  for  the  simulations.  It  also  sends  a  “step-list”  to  the 
CAM-8  to  specify  the  sequence  of  kicks  and  scans  that  evolve  the  lattice- 
gas  in  time.  One  can  view  the  lattice-gas  simulation  in  real-time  since  a 
custom  video  module  captures  site  data  for  display  on  a  VGA  monitor,  a 
useful  feature  for  lattice-gas  algorithm  development,  test  and  evaluation.  The 
CAM-8  has  built-in  25-bit  event  counters  allowing  real-time  measurements 
without  slowing  the  lattice-gas  evolution.  This  feature  is  used  to  do  real¬ 
time  coarse-grain  block  averaging  ol  the  lattice-gas  number  variables  and 
to  compute  the  components  of  the  momentum  vectors  for  each  block.  The 
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amount  of  coarse-grained  data  is  sufficiently  small  to  be  transferred  back  to 
the  front-end  host  for  graphical  display  as  an  evolving  flow  field  within  an 
X- window. 

3  Lattice-Gas  Automaton:  An  Exactly  Com¬ 
putable  Dynamical  System 

A  Boolean  formulation  of  an  exactly  computable  dynamical  system,  known 
as  a  lattice-gas,  may  be  stated  in  a  way  that  is  consistent  with  the  Boltzmann 
equation  for  kinetic  transport.  In  essence  the  lattice-gas  dynamics  are  a  sim¬ 
plified  form  of  molecular  transport  as  we  restrict  ourselves  to  a  cellular  phase 
space.  The  macroscopic  equations,  in  particular  the  continuity  equation  and 
the  Navier-Stokes  equation,  are  obtained  by  coarse-graining  over  a  discrete 
microdynamical  transport  equation  for  a  number  of  Boolean  variables.  The 
scheme  employs  the  finite-point  group  symmetry  of  a  crystallographic  spatial 
lattice.  It  is  somewhat  inevitable  that  to  obtain  an  exactly  computable  rep¬ 
resentation  of  fluid  dynamics  one  must  perform  a  statistical  treatment  over 
discrete  number  variables. 

Before  the  basic  lattice-gas  microdynamical  transport  equation  is  intro¬ 
duced,  some  notational  conventions  are  needed.  Consider  a  spatial  lattice 
with  N  total  sites.  The  fundamental  unit  of  length  is  the  size  of  a  lattice 
cell,  /,  and  the  fundamental  unit  of  time,  r,  is  the  time  it  takes  for  a  speed-one 
particle  to  go  from  one  lattice  site  to  a  nearest  neighboring  site.  Particles, 
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•wath  unit  mass  m,  propagate  on  the  lattice.  The  unit  lattice  propagation 
speed,  is  denoted  by  c  —  Particles  occupy  this  discrete  space  and  can 
have  onh'^  a  finite  number  B  of  possible  momenta.  The  lattice  vectors  are 
denoted  by  where  a  =  l,2,...,B.  For  example,  for  a  single-speed  gas  on 
a  triangular  lattice,  o  =  1, 2, . . . ,  6.  A  particle’s  state  is  completely  specified 
at  some  time,  t,  by  specifying  its  position  on  the  lattice,  Xi,  and  its  momen¬ 
tum,  Pi  =  mceai,  at  that  position.  The  particles  obey  Pauli  exclusion  since 
only  one  particle  can  occupy  a  single  momentum  state  at  a  time.  The  total 
number  of  configurations  per  site  is  2^.  The  total  number  of  possible  single 
particle  momentum  states  available  in  the  system  is  A'toiai  =  BA' .  W'ith  P 
particles  in  the  system,  denote  the  filling  fraction  by  d  = 

The  number  variable,  denoted  by  na(x,t),  takes  the  value  of  one  if  a 
particle  exists  at  site  x  at  time  t  in  momentum  state  me  Ca,  and  takes  the 
value  of  zero  otherwise.  The  evolution  of  the  lattice-gas  can  then  be  written 
in  terms  of  Uq  as  a  two-part  process:  a  collision  part  and  a  streaming  part. 
The  collision  part  reorders  the  particles  locally  at  each  site. 

n'^(x,t)  =  na(x,t)  -l-Oa(n(x,t)),  (1) 

where  fla  represents  the  collision  operator  and  in  general  depends  on  all  the 
particles,  n  at  the  site.  So  as  a  shorthand  suppress  the  index  on  the  occupa¬ 
tion  variable  when  it  is  an  argument  of  fia[n(x,t)]  to  represent  this  general 
dependence.  In  the  streaming  part  of  the  evolution  the  particle  at  position 
X  “hops”  to  its  neighboring  site  at  x  -I-  /ea  and  then  time  is  incremented  bv 
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(2) 


n|^(x  +  /ea,t +  r)  =  na(x,t)  +  na[n(x,t)]. 

Equation  (2)  is  the  lattice-gas  microdynamical  transport  equation  of  motion. 

The  collision  operator  can  only  permute  the  particles  locally  on  the  site 
since  we  wish  the  local  particle  number  to  be  conserved  before  and  after  the 
collision.  That  is, 

n(x,t)  =  ^n;(x,t)  =  X^na(x,t).  (3) 

a  CL 

Equation  (3)  defines  the  local  number  density.  Summing  Eq.  (2)  over  lattice 
directions  then  implies  the  constraint  on  the  collision  operator, 

a 

We  may  define  the  local  momentum  as 

Pi(x,t)  =  mc^eain'^{^,t)  =  mc52eat^a(x,t),  (5) 

a  0. 

which  of  course  must  also  be  conserved  before  and  after  a  collision.  Again, 
this  imposes  a  constraint  on  the  collision  operator. 

=  (6) 

a 

As  a  matter  of  notation  it  should  be  understood  that  whenever  a  directionally 
dependent  quantity  is  written,  its  subscripted  index  is  taken  modulo  B.  Using 
the  number  variable,  for  example,  it  is  understood  that 

TT'a+b  ~  ^modB(a+6)-  ('^) 
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As  a  shorthand,  a  negative  index  will  represent  the  antiparallel  direction,  so 
since  =  -e<j  we  may  write 

=  (S) 

4  Isotropic  Lattice  Tensors 

\Ye  construct  an  rank  tensor  composed  of  a  product  of  lattice  vectors  [7] 

(9) 

a 

where  a  =  1,. . .  ,B.  All  odd  rank  E  vanish.  We  wish  to  express  in 

terms  of  Kronecker  deltas,  %  =  1  for  i  =  j  and  zero  otherwise.  We  can  turn 
this  problem  of  expressing  the  i^-tensors  in  terms  of  products  of  Kronecker 
deltas  into  a  problem  of  combinatoric  counting.  W'e  use  the  following  tensors: 

Aj  =  k  (10) 

^tjkl  —  ^ij^k!  +  ^ik^jl  +  ^il^kj  (^1) 

and  so  forth.  Then  we  know  that  if  E  is  isotropic  it  must  be  proportional  to 
A,  therefore 

a  (12) 

and  that  the  constant  of  proportionality  may  be  obtained  by  counting  the 
number  of  ways  we  could  write  a  term  comprising  a  product  of  n  Knonecker 
deltas.  Consider  for  example  the  case  n  =  2.  Since  the  Knonecker  delta  is 
symmetric  in  its  indices,  the  following  four  products  are  identical:  5ij5ki  = 
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5ij5ik  =  SjiSki  =  6ji6ik.  The  degeneracy  is  2^.  Furthermore,  the  order  of  the 
Kronecker  deltas  also  doesn’t  matter  since  they  commute;  that  is,  = 

5ki5ij.  The  degeneracy  is  2!.  For  the  case  where  n  is  arbitrary,  there  are 
2"  identical  ways  of  writing  the  product  of  n  Kronecker  deltas.  For  each 
choice  of  indices,  there  are  an  additional  n!  number  of  ways  of  ordering  the 
products.  Therefore,  the  total  number  of  degeneracies  equals  2"n!  =  (2n)!!. 
The  total  number  of  permutations  for  2n  indices  equals  (2n)!.  So  from  this 
counting  procedure  we  know  that  consists  of  a  sum  of  ^2^)” 
terms. 

The  following  relations  will  be  very  useful  throughout  later  developments 


E^ 


0 

0 


E^ 


B 

D{D  +  2) 


{5ij5ki  +  SikSji  +  SiiSkj) 


In  general,  the  lattice  tensors  are 


E^n+l  ^  Q 

p2n  _  _ B _ ^2n 

D{D  +  2)-- -{0  +  271-2) 


(13) 

(14) 

(15) 

(16) 


(17) 

(18) 
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5  Triangular  Lattice 

In  a  triangular  lattice  there  are  six  vectors,  enumerated  in  this  section  b}-  the 
convention 

©a  =  (cos  Y,-sin  y)  ,  (19) 

where  a  =  1, 2, . . . ,  6.  The  spatial  coordinates  of  the  lattice  sites  may  be 


(1,5)  *r - (2,5)t - (3,5)*: - (‘•.SK; - (5,5)^| - j 

.  (1,4)V-— (2,4)V-  — (3,4]>— -(Mii--— (5,4jv 

■  1  ^  ^  /  *  /  /  / 


\ 


y 


a  =  6  a  =  3  (1.3)V--(2,3)(---(3,3)^-(4,3)  Y--(5.3K^ 

,  5"  v.-""  a  =  4  j  (1 ,2) y  -  —  (2^y-  -  -  (3,2)>-  - «  -  {4,2)y  -  -  -  (5^> 


(a) 


(b) 


Figure  3:  Triangular  Lattice  Convention:  (a)  Lattice  vector  label  convention; 
(b)  Triangular  lattice  convention  with  lattice  directions  a  =  3  up  and  a  =  6 
down.  Coordinates  above  the  lattice  nodes  are  {i,j)  memory  array  indices. 

expressed  as  follows 

Xy  =  (*  “  mod  2),  (20) 

where  i  and  j  are  rectilinear  indices  that  specify  the  data  memory  array 
location  used  to  store  the  lattice-gas  site  data. 
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Let  s  =  (7  mod  2) (r  mod  2) .  Given  a  particle  at  site  (lj),  it  mav  be 
shifted  to  a  site  r  lattice  units  away  to  a  remote  site  {i'J')  by  the  mapping 


(‘+  2 

(21) 

= 

{1-5 r) 

(22) 

(*■',. ?")3  = 

(23) 

= 

(i- 5 

(24) 

= 

(25) 

= 

{i  +  r,j) 

(26) 

where  denotes  the  shifted  site,  that  is,  {i,j)  ^  with  a  shift 

vector  T  —  tGq  and  where  division  by  2  is  considered  integer  di\'ision. 

6  Local  Collision  Rules 


In  two  dimensions  we  may  use  a  triangular  lattice,  with  six  bits  per  site 
encoding  the  occupation  numbers  of  the  six  possible  momentum  states.  Let 
Ua  be  the  input  bits  and  n'„  be  the  output  bits  of  a  local  collision.  A  general 
collision  operator  is  constructed  as  follows 


(27) 

{Cd 

where  {G}  is  a  set  of  occupied  particle  states,  a  =  ±1  is  a  scalar  coefficient, 
and  where  each  term  in  the  sum  is  written  in  factorized  form  as 


QaiHi  •  ■  ■  1  ^k) 


B 


l-na 


-fil 


1  /i'a+U-  j=l 


(28) 
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Table  1:  Simple  Right-Handed  Collision  Table 


no 

ni 

n2 

n3 

n4 

no 

n(, 

n\ 

n'2 

n4 

^5 

1 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

1 

0 

1 

0 

1 

0 

Then  the  FHP  collision  operator  is  the  following: 

nfP  =  ^(5a(l,4)  -h  iQa(2,5)  -  Q„(0,3)  +  Qa(l,3,5)  -  Q„(0,2,4)  (29) 

or  for  a  =  0  this  is  explicitly 

Qo  =  ^nin4(l  —  no)(l  —  n2)(l  —  n3)(l  —  ns)  + 

^n2n5(l  -  n())(l  -  ni)(l  -  n3)(l  -  n^)  - 
non3(l  -  ni)(l  -  n2)(l  —  n4)(l  -  ns)  -t- 
nin3n5(l  -  no)(l  -  n2)(l  -  n4)  - 
non2n4(l  -  ni)(l  -  n3)(l  -  n^). 

Note  that  it  is  sufficient  to  give  only  fio  since  the  other  components  of  the 
collision  operator  can  be  obtained  simply  by  incrementing  the  indices  of  flo 
owing  to  the  six-fold  symmetry  of  the  collisions.  The  factors  of  one-half  in 
Eq.  (29)  are  transition  probabilities  for  the  2-body  collisions,  indicating  a 
coin  toss  is  performed  to  choose  betrveen  even  or  odd  chirality. 
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Table  2;  Simple  Left-Handed  Collision  Table 


no  ni  n2  Ui  ns 

n'o  n[  n’2  n'2  n'^  n'g 

10  0  10  0 
010  0  10 

0  0  1  0  0  1 

0  10  0  10 

0  0  1  0  0  1 
10  0  10  0 

10  10  10 

0  10  10  1 

0  10  10  1 
10  10  10 

100100 ->001001  I 

01  001  0->1  001  00 

1  001001 ->010010 

Even  Chirality 

1  001  00->  01  001  0 

i  01001  o->  001  001  1 

i  \  0  0  1  0  0  1  ->  1  0  0  1  0  0 

i  +  ^ 

Odd  Chirality 

1  101010 ->010101 

T 

01  01  01 ->1  01  01  0 

Figure  4:  FHP  Collision  Rules.  Enumeration  of  FHP  2-body  collisions,  even 
and  odd  chirality,  and  3-body  collisions. 

The  possible  two-body  and  three-body  collisions  represented  by  Eq.  (29) 
are  illustrated  in  Figure  4.  For  two-dimensional  flow,  there  are  four  invari¬ 
ants,  the  mass,  two  components  of  the  momentum,  and  the  energy.  With 
only  the  2-body  collision  in  Figure  4,  there  is  an  additional  invariant;  the 
difference  in  the  particle  number  along  each  of  the  three  lattice  directions. 
The  3-body  collisions  in  Figure  4  were  include  in  the  FHP-model  to  remove 
this  spurious  invariant.  Consequently,  the  collisions  enumerated  in  Figure  4 
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are  the  minimally  sufficient  set  to  produce  hydrodynamic  behavior  in  the 
continuum  limit. 

7  Long-Range  2-Body  Interactions 

An  interparticle  potential,  V^(x  —  x^),  acts  on  particles  spatiallv  separated 
by  a  fixed  distance,  x  -  x'  =  r.  An  effective  interparticle  force  is  caused 
by  a  non-local  exchange  of  momentum.  Momentum  conservation  is  violated 
locally,  yet  it  is  exactly  conserved  in  the  global  dynamics. 


.....  /  .  \ 

x’  x\  /x’ 


(b)  (c) 

Figure  5:  Bounce-Back  and  Clockwise  Bound  States.  Simple  bound-state 
orbits  due  to  a  long-range  attractive  interaction  where  the  dotted  path  in¬ 
dicates  the  particle’s  closed  trajectory:  (a)  partition  directions:  (b)  bounce- 
back  orbit  with  |Ap|  =  2  and  zero  angular  momentum;  and  (c)  clockwdse 
orbit  with  |  Ap|  =  1  and  one  unit  of  angular  momentum.  Head  of  the  dashed 
arrows  indicates  particles  entering  the  sites  of  partition  tq  at  time  t.  Tail  of 
the  black  arrow^s  indicates  particles  leaving  those  sites  at  time  t  -t-  r.  The 
counter-clockwise  orbit  is  not  shown. 


(a) 


For  an  attractive  interaction,  there  exists  bound  states  in  which  two  par¬ 
ticles  orbit  one  another.  Since  the  particle  dynamics  are  constrained  by  a 
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Table  3:  Lattice  Vector  Components 


a  x-component 

y-component 

0  -1 

0 

1  -i 

2  - 

2 

3  1 

0 

4  “ 

^  2 

x/3 

2^ 

r  1 

x/3 

^  2  

2 

crystallographic  lattice  we  expect  polygonal  orbits.  In  Figure  5  we  have 
depicted  two  such  orbits  for  a  hexagonal  lattice-gas.  The  range  of  the  in¬ 
teraction  is  r.  Two-body  single  range  attractive  interactions  are  depicted 
in  Figures  5b  and  5c,  the  bounce-back  and  clockwise  orbits  respectively. 
Momentum  exchanges  occur  along  the  principal  directions.  The  interaction 
potential  is  not  spherically  symmetric,  but  has  an  angular  anisotropy.  In 
general,  it  acts  only  on  a  finite  number  of  points  on  a  shell  of  radius  |.  The 
number  of  lattice  partitions  necessary  per  site  is  half  the  lattice  coordina¬ 
tion  number,  since  two  particles  lie  on  a  line.  Though  microscopically  the 
potential  is  anisotropic,  in  the  continuum  limit  obtained  after  coarse-grain 
averaging,  numerical  simulation  indicates  isotropy  is  recovered. 

8  A  Simple  Example:  Bounce-Back  Orbit 

A  long-range  lattice-gas  of  the  type  we  are  considering  still  possesses  the 
usual  local  dynamics  of  a  hydrodynamic  lattice-gas.  To  extend  the  local 
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lattice-gas  update  rules  to  include  long-range  interactions,  we  use  two  addi¬ 
tional  bits  of  local  site  data.  This  will  allow  us  to  implement  a  long-range 
interaction  using  strictly  local  updating  and  therefore  the  algorithm  remains 
'’embarrassingly”  parallel  just  as  a  usual  local  lattice-gas.  The  two  addi¬ 
tional  bits  will  denote  the  occupation  numbers  of  messenger  particles,  or 
“photons” .  The  idea  of  using  messenger  particles  was  introduced  by  Appert 
et  a/. [8].  \^'e  have  two  types  of  messenger  states,  to  represent  incoming  and 
outgoing  conditions,  and  we  denote  the  messengers  as  zi  and  Zr- 

For  the  simplest  long-range  lattice-gas  model,  we  therefore  use  eight  bits 
of  local  site  data.  Since  long-range  interactions  occur  between  remote  spa¬ 
tial  sites,  say  x  and  f',  the  messenger  particles  will  travel  either  parallel  or 
antiparallel  to  the  vector  f=x-x\  All  pairs  of  sites  throughout  the  entire 
space  that  are  separated  by  the  vector  f  can  therefore  all  be  updated  in  paral¬ 
lel.  We  refer  an  update  step  of  all  pairs  of  2-body  interactions  along  direction 
r  as  a  partition,  denoted  by  Fr-  All  possible  two-bod]y  interaction  pairs  are 
then  computed  by  performing  all  possible  partitions  of  the  space.  For  this 
reason,  it  requires  many  scans  for  the  space  to  perform  a  single  long-range 
interaction  step. 

In  our  two-dimensional  example  using  a  triangular  lattice,  there  are  three 
possible  partitions.  The  number  of  partitions  is  ne’ver  smaller  than  half  the 
lattice  coordination  number.  In  the  two-dimensional  case,  the  simplest  long- 
range  lattice-gas  algorithm,  though  perhaps  not  the  most  efficient  algorithm, 
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is  to  use  three  sequential  scans  of  the  space,  each  scan  performing  the  up¬ 
dating  necessary  for  a  single  partition.  See  Figure  5a.  Often,  depending  on 
the  complexity  of  the  long-range  interactions  and  the  dimensionality  of  the 
lattice,  it  is  possible  to  perform  simultaneous  updating  of  multiple  partitions. 
This  of  course  is  more  desirable,  yet  it  causes  more  complexity.  Furthermore, 
this  updating  requires  an  extra  pair  of  messenger  particles  for  each  partition 
to  be  simultaneously  updated.  For  simplicity,  we  will  not  deal  with  this  case 
here,  however  our  implementation  on  the  CAM-8  does  use  simultaneous  par¬ 
tition  updating— repulsive  and  attractive  partitions  are  performed  in  parallel 
using  four  messenger  bits. 

Let  us  consider  a  simple  example  of  the  long-range  lattice-gas  algorithm, 
the  minimal  model  of  Appert.  Here  we  consider  only  bounce-back  attractive 
interactions.  Suppose  there  is  a  single  particle  at  site  f  =  0  and  there  is  also 
a  single  particle  at  site  x'  =  ri:  that  is,  no{x)  =  1,  ns{x)  =  0,  n.o(x )  =  0 
and  nsix')  =  1  with  all  other  bits  at  x  and  x'  being  zero.  See  Figure  5b  for 
a  diagram  of  this  situation.  Here  we  are  using  the  bit  convention  shown  in 
table  3.  Then  the  two  particles  are  separated  by  a  distance  r  and  are  moving 
away  from  each  other.  The  attractive  long-range  interaction  will  effectively 
flip  their  respective  directions  making  no(f)  =  0,  nz{x)  =  1,  no(^)  =  1  and 
n3(x')  =  0  so  that  the  two  particles  will  now  be  moving  toward  each  other. 
There  is  a  local  momentum  change  of  2mci  at  and  an  opposite  momentum 
change  of  — 2mci  at  x.  Locally  momentum  is  not  conserved,  but  nonlocally 
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it  is. 


The  first  step  of  the  long-range  interaction  sequence  is  to  choose  a  parti¬ 
tion,  say  Fr,  and  then  to  emit  messenger  particles  along  the  partition  axis. 
The  basic  local  rule  for  this  first  step  is  the  following:  a  photon  is  emitted  at 
a  site  if  there  exists  a  particle  at  that  site  that  can  partake  in  a  long-range 
interaction.  Another  way  of  expressing  this  rule  is:  send  only  if  you  can 
Tccezve.  Obviously,  for  a  particle  to  partake  in  an  interaction  there  must  be 
both  a  particle  and  a  hole  at  that  site.  The  factorized  probability  of  having 
such  a  situation  is  just  d(l  —  d).  So  to  continue  with  our  example,  for  a 
photon  to  be  emitted  at  some  site  x  parallel  or  antiparallel  to  a  partition 
direction  i,  we  use  the  following  rule 

Zr{x)  =  no{x){l  -  n^{x))  (30) 

zi{x)  =  n3(f)(l  -  no(f)).  (31) 

Note  that  according  to  this  local  rule,  only  one  photon  can  be  created  at  a 
site,  and  consequently  w^e  eliminate  the  possibility  of  a  long-range  interaction, 
say  of  range  2r,  mediated  through  a  doubly  occupied  site.  For  two  sites 
separated  by  the  interaction  distance  r,  the  important  consequence  of  the 
emission  step  is  that  if  both  sites  send  photons,  both  wall  necessarily  receive 
them,  which  strictly  enforces  nonlocal  momentum  conservation.  Give  and  ye 
shall  receive  (provided  yours  is  received).  Letting  Za  =  and  Z-a  =  zi/m 
general  w-e  can  write  the  emission  step  of  the  minimal  interaction  as 

z^{x)  =  n_„(f)(l  -  n„(f)),  (32) 
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where  a  =  0, 1,2  covers  all  the  partitions. 

After  the  emission  step,  a  long-range  kick  of  the  messenger  bits  follows. 
In  the  simple  example,  all  photons  zi  are  kicked  along  -ri  and  all  photons 
.2,.  are  kicked  along  ri.  In  general  for  the  long-range  kick  we  have 

z'^{x  +  rea)  =  Za{x).  (33) 

Finally,  we  have  the  absorption  step  of  the  long-range  interaction  sequence. 
Here  the  local  particle  momentum  state  is  updated  as  the  particles  flip  their 
directions  in  our  example 

n'^{x)  =  -F  Zi{x)nQ{x){\  —  Ti2i{x))  —  Zj.{x)n^{x)(l  —  no^x)) 

(34) 

=  uq^x)  —  tiq{x))  —  Zi[x)nQ{x){l  n^{x)) . 

(35) 

Aloreover,  in  this  step  all  the  messenger  bits  are  set  to  zero  throughout  the 
entire  space.  For  any  direction,  the  local  absorption  rule  could  be  written 
more  simply  as 

n(,(f)  =  na(f)  +  z'_^{x)za{x)  -  z'^{x)z^a{S).  (36) 

Substituting  in  Eqs.  (32)  and  (33)  into  Eq.  (36),  we  have  a  single  Boolean 
expression  in  terms  of  number  variables  for  a  single  long-range  interaction 
step  for  partition  F^  as  follows 

n'^{x)  =  na{x)  -f 
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Table  4:  Long-range  Interaction  Sequence 


events 

Uq  (^) 

Zi{x) 

Zf.  (^) 

Tla{^  ) 

Zl{x') 

Zr  ) 

initial 

100000 

0 

0 

000100 

0 

0 

emit 

100000 

0 

1 

000100 

1 

0 

kick 

100000 

1 

0 

000100 

0 

1 

absorb 

000100 

0 

0 

100000 

0 

0 

na{x  +  rea)(l  —  n^ai^  +  ?’ea))^-a(^){l  ^a(^)) 

n_a(x  -  rea)(l  -  ne(f  -  re„))n<,(f)(l  -  n^a{x)) 

(37) 

For  convenience  we  define  a  long-range  collision  operator,  Pa,  as  follo-ws 


Pa(f)  =  z'_^{x)Za{x), 


(38) 


SO  that  we  may  write 

n(^(f)  =  na{x)  -i-  Pa{x)  -  P-a{x).  (39) 

The  state  data  for  this  simple  example  we  have  been  considering  are  given 
in  Table  4,  which  represents  all  the  steps  of  a  long-range  interaction  sequence 
for  a  partition  along  the  x-axis. 

9  Another  Example:  Clockwise  Orbit 

To  continue  illustrating  our  implementation  of  a  long-range  lattice-gas,  in 
this  section  we  again  consider  a  system  with  a  single  attractive  interaction  of 
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range  r,  however  the  local  momentum  states  participating  in  the  interaction 
are  not  along  the  partition  direction.  However,  in  the  example  given  here,  the 
momentum  exchange  is  still  along  the  partition  direction  so  that  the  long- 
range  interaction  remains  a  central-body  one,  resulting  in  a  bound  state 
with  two  particles  trapped  in  a  clockwise  orbit.  (Note  that  the  restriction  to 
central-body  forces  is  not  necessary,  but  is  presented  here  for  convenience.) 
In  this  slightly  more  complicated  example,  the  local  rules  for  photon  emission 
and  absorption,  Eqs.  (32)  and  (36)  respectively,  have  a  more  general  form 
with  the  implication  that  the  emission  and  absorption  of  photons  is  different 
from  the  previous  example  of  the  bounce-back  orbit,  and  the  difference  should 
be  noted  when  making  look-up  tables  to  do  this  computation.  The  local 
photon  emission  rules  can  be  written 

Za{x)  =  nc(f)(l  -  n(i(f))  (40) 

z-a{x)  =  ng{x){\  -  nh{x))  (41) 

where  the  bits  c,  d,  .g,  and  h  must  by  chosen  so  momentum  is  conserved 

ec  -  -h  Cg  -  e/,  =  0  (42) 

as  well  as  be  constrained  by  central-body  parallel  and  perpendicular  momen¬ 
tum  exchange  conditions 

(ec  -ed-  eg  +  h)  ■  r  -  2Ap 

{ic  —  —  €g  +  ih)  X  f  —  0, 
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(43) 

(44) 


where  Ap  is  the  momentum  change  per  site  due  to  the  long-range  interaction. 
Equations  (40)  and  (41)  differ  for  Eqs.  (30)  and  (31)  for  the  bounce-back 
orbit  by  allowing  two  photons  to  be  emitted  at  a  single  site. 

To  be  explicit,  for  the  two-dimensional  triangular  lattice,  we  can  satisfy 
Eqs.  (42),  (43),  and  (44)  by  choosing  the  indices  c,d,gji  as  follows: 


c  =  a-2  (45) 

d  =  a- 1  (46) 

p  =  -c  (47) 

h  =  -d.  (48) 

An  example  of  this  choice  of  indices  is  illustrated  in  Figure  5c.  Then  the 
emission  rules,  Eqs.  (40)  and  (41),  are  simply 

Za{.x)  =  na-2{x)[l  “  na-i{x)]  (49) 


Since  the  kicking  of  the  photons  is  the  same  in  this  example  as  in  the  previous 
one,  Eq.  (33)  still  holds 

4(f  •+•  rCa)  =  Zaix). 

By  re-expressing  Eq.  (36)  more  generally,  we  can  wuite  a  local  absorption 
rule 

n'ai^)  =  Mx)  +  z'_^^^^){x)Za+l{x)  -  (f  )z_(a_l)  (^)  (50) 

or  more  elegantly 

n(,(f)  =  na{x)  +  Pa+l{x)  -  P-a+l{x).  (51) 
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Table  5;  Long-Range  Interaction  Sequence  with  Two  Photons  Emitted  at  a 
Single  Site 


events 

na{x) 

Zi{x) 

Zr  (^) 

) 

Zi{x') 

Zr  ) 

initial 

010010 

0 

0 

000010 

0 

0 

emit 

010010 

1 

1 

000010 

1 

0 

kick 

010010 

1 

0 

000010 

0 

1 

absorb 

001010 

0 

0 

000001 

0 

0 

Substituting  from  Eqs.  (49)  and  (33)  into  Eq.  (50)  and  after  some  manipu¬ 
lation  of  the  indices,  we  have  a  single  boolean  expression  in  terms  of  number 
variables  for  a  single  long-range  interaction  step  for  partition  Pr  as  follows 

n(j(f)  =  na{x)  -b 

na+2{x  +  rea+i)[l  -  n_a(f  +  rea+i])n„_i(x)[l  -  na{x)\  - 
n_a(f  -  re„_i)[l  -  Ua-^ix  -  rea-i)]na{x)[l  -  n„+i(f)]. 

Table  5  gives  the  local  site  data  for  the  x-axis  partition  of  a  clockwise 
orbit.  The  particle  724  (f)  acts  as  a  kind  of  spectator  is  this  example,  illus¬ 
trating  that  two  photons  can  be  emitted  from  a  single  site.  It  is  also  possible 
to  have  tw'o  photons  absorbed  at  a  single  site.  Let  us  consider  a  back-to- 
back  interaction  over  three  sites.  Suppose  there  are  particles  at  sites  x  =  0, 
x'  =  ri,  and  x"  =  2ri.  Table  6  gives  the  site  data  for  these  sites  where  there 
are  two  photons  emitted  and  absorbed  at  x'  in  the  middle  location. 

The  minimal  model,  using  only  an  attractive  interaction,  models  a  fluid 
with  liquid  and  gas  phases  at  zero  temperature.  Figure  6  shows  the  time 
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Table  6:  Long-Range  Interaction  Sequence  with  Two  Photons  Emitted  and 
Absorbed  at  Site  x'  in  a  Back-to-Back  Interaction 


events 

na{x) 

zi{x) 

Zr{x) 

na{x') 

Zi{x') 

Zr{x') 

na{x'') 

Zl{x'') 

Zr{x") 

initial 

010000 

0 

0 

010010 

0 

0 

000010 

0 

0 

emit 

010000 

0 

1 

010010 

1 

1 

000010 

1 

0 

kick 

010000 

1 

0 

010010 

1 

1 

000010 

0 

1 

absorb 

001000 

0 

0 

001001 

0 

0 

000001 

0 

0 

evolution  of  the  phase  separation  process  in  this  case  at  a  density  d  =  0.07 
and  interaction  range  r  =  Ql  and  illustrates  the  type  of  physical  simulation 
that  can  be  achieved  with  the  simplest  long-range  lattice-gas  algorithm. 
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Figure  6:  [Liquid-Gas  Phase  Separation.  Time  evolution  of  a  liquid-gas  phase 
separation  for  a  lattice-gas  with  long  range  attractive  interactions  at  range 
r  =  6  on  a  1024  x  1024  lattice  starting  with  a  uniformly  random  configuration 
of  density  d  —  0.07. 


27 


10  Symmetries  in  Long  Range  Interactions 

For  an  attractive  interaction  as  mentioned  earlier,  we  would  expect  that  there 
exists  a  bound  state  in  which  two  particles  orbit  one  another,  for  example 
the  clockwise  orbit  shown  in  Figure  5c.  The  range  of  the  interaction  is  r. 
This  hexagonal  orbit  is  also  shown  in  the  bottom  right  corner  of  Figure  7. 
In  this  figure  the  hexagon’s  radius  is  also  labeled  as  r  and  should  not  be 
confused  with  the  interaction  range,  that  is  the  hexagon  s  diameter.  Simlarlv 
a  counter-clockwise  orbit  is  possible.  See  the  top  left  corner  of  Figure  7. 
A  time-reversal  invariance  exists  between  these  two  cases  with  respect  to 
conjugation  and  parity  operations,  as  depicted  in  Figure  7. 

These  diagrams  describe  the  possible  2-body  collisions  that  can  be  so 
generated,  including  repulsive  interactions.  This  logical  correspondence  be¬ 
tween  the  different  types  of  2-body  interactions  allows  one  to  achieve  a  more 
efficient  implementation  of  a  long-range  lattice-gas  algorithm  than  what  we 
have  achieved  on  the  CAM-8,  since  we  have  not  used  this  form  of  logical  econ¬ 
omy.  Furthermore,  because  there  is  a  correspondence  property  for  r  =  0,  the 
situation  reduces  to  the  2-body  collisions  in  the  FHP  lattice-gas 

The  long  range  interactions  considered  here  have  the  following  properties 
that  simplify  a  computational  implementation:  1)  there  exists  onlv  parallel 
momentum  exchange,  a  restriction  for  modelling  central  body  forces;  2)  the 

^For  r  =  0,  the  |Ap|  =  1  interactions  reduce  to  a  rotation  of  the  states,  and 

and  the  |Apl  =  2  interactions  reduce  to  the  identity  operation. 
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e  <-  iB 


Figure  7:  Long-Range  2-Body  Interaction  Terms.  Examples  of  two-body 
finite  impact  parameter  collisions  along  the  ro-direction.  The  four  terms  of 
the  interaction  Hamiltonian  are  for  |Ap|  =  1.  Input  states  are  depicted  by 
dashed  arrows  and  output  states  are  by  black  arrows. 
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interaction  acts  only  along  the  principal  lattice  directions;  3)  there  exists 
time-reversal  invariance  with  respect  conjugation  and  parity  operations;  and 
4)  a  bias  to  the  interactions  can  be  assigned  by  coupling  to  a  heat-bath 
reservoir.  The  previous  examples,  bounce-back  and  clockwise  orbits,  has 
been  discussed  properties  1  and  2.  Property  3  has  been  discussed  in  this 
section.  Property  4  has  been  discussed  elsewhere  [9]. 

11  Central-Body  Interaction  Neighborhood 
and  2D  Crystallization 

In  the  previous  two  sections,  two  examples  of  the  long-range  interactions, 
the  bounce-back  and  clockwise  orbits  on  a  two  dimensional  triangular  lattice 
were  illustrated  .  In  general,  the  long-range  interaction  step  will  involve  many 
partitions,  both  attractive  and  repulsive  interactions,  and  multiple  ranges. 
In  my  CAM-8  implementation  of  a  long-range  lattice-gas  with  central-body 
interactions,  I  use  12  neighbors  in  two  dimensions,  as  indicated  in  Figure  8. 
The  triangular  lattice  is  superposed  over  a  square  lattice,  which  appears 
rhomboidal  in  the  figure.  The  square  lattice  is  often  used  for  embedding  the 
site  data  into  computer  memory,  which  is  rectilinear.  This  kind  of  embed¬ 
ding  is  the  simplest  and  used  for  simulations  that  possess  periodic  boundary- 
conditions.  The  reason  for  using  12  neighbors  is  to  try  to  achieve  a  higher 
degree  of  local  s}^mmetry.  In  doing  molecular  d^'^namics  modeling  with  a 
multi-long-range  lattice-gas,  we  have  found  that  12  neighbors  are  necessarj'^ 
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(-R,-2R) 


Figure  8:  Twelve  Neighbors  on  a  Triangular  Lattice  that  can  Participate  in 
Long-Range  Central-Body  Interactions.  Interaction  range  D  =  7  is  depicted, 


where  i?  =  ~  4  to  within  1.03  percent  error. 


Computer  memory  space 


coordinates  (x,  y)  are  given  adjacent  to  each  neighboring  site. 


to  recover  macroscopic  isotropy.  In  particular,  12  neighbors  are  necessary  to 
have  the  emergent  crystalline  solid  be  able  to  freely  rotate  in  space.  A  mean- 
field  analysis  of  the  lattice-gas  crystallization  method  has  been  presented 
elsewhere  [10].  Figure  8  shows  a  ring  of  range  7  lattice  spacings. 

To  implement  the  crystallization  algorithm,  we  use  up  to  eight  ranges 
in  two  dimensions,  that  is,  eight  rings  of  the  type  shown  in  Figure  8  for 
a  total  of  96  neighbors.  Half  the  rings  are  used  for  attractive  interactions 
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t=500 


t=5000 


Hexagonal 
/Close  Pack 


512  X  512  Lattice 
6  ranges 
12  directions 


Figure  9:  Time  Evolution  of  Crystallization  in  a  Two  Dimensional  Lattice- 
Gas  with  Multiple  Fixed-Range  2-body  Interactions.  The  resulting  crystal  is 
in  a  hexagonal-close- pack  configuration  since  the  coarse-grained  interatomic 
potential  is  radially  symmetric.  The  underlying  lattice  is  512  x  512.  Started 
with  a  uniformly  random  configuration  at  d  =  0.1.  Twelve  directions  are 
used  for  long-range  momentum  exchanges.  Grain  boundaries  and  defects  are 
observed  during  the  early  stages  of  the  crystal  formation. 


and  the  other  half  are  used  for  repulsive  interactions.  Typically,  the  inner 
rings  are  attractive,  the  middle  rings  are  repulsive,  and  the  outer  rings  are 
again  attractive.  Since  four  photon  bits  are  used  in  our  implementation,  and 
since  each  ring  is  either  attractive  or  repulsive,  two  rings  are  affected  by  a 
simultaneous  partitioning  of  the  space.  For  the  attractive  interaction,  there 
are  five  types  of  orbits:  bounce-back  with  Ap  =  2,  clockwise  and  counter¬ 
clockwise  with  Ap  =  1,  and  clockwise  and  counter-clockwise  with  Ap  =  -\/3. 
Consequently,  there  are  five  types  of  repulsive  interactions,  which  are  just 
the  conjugates  of  the  five  attractive  ones.  Since  there  are  three  partition 
directions  for  a  triangular  lattice,  it  takes  5*3=15  partitions  of  the  space 
to  completely  update  an  attractive  ring  and  a  repulsive  ring  simultaneously. 
To  compute  8  rings  therefore  takes  4*15=60  scans  of  the  space.  Therefore, 
since  the  local  collisions  require  a  single  scan,  it  takes  a  total  of  61  scans  to 
complete  one  time  step. 

A  tw-o  dimensional  example  using  six  interaction  ranges,  Avith  an  un¬ 
derlying  512  X  512  lattice,  of  this  time-dependent  crystallization  process  is 
given  in  Figure  9  and  illustrates  the  type  of  molecular  dynamics  simulation 
that  can  be  achieved  Avith  a  more  complex  long-range  lattice-gas  algorithm. 
The  resulting  crystal  is  in  a  hexagonal-close-pack  configuration  since  we  haA^e 
striA'-ed  to  make  the  coarse-grained  interatomic  potential  be  radially  symmet¬ 
ric.  This  long-range  lattice-gas  model  had  six  interaction  ranges:  r  =  -2 
,-7,  19,  21,  —24,  -26.  Here  the  negatiA’e  sign  preceding  the  range  denotes 


33 


an  attractive  interaction  at  that  range. 


12  Conclusion 

It  should  be  noted  that  although  the  lattice-gas  molecular  dynamics  algo¬ 
rithm  that  we  have  described  above  requires  61  scans,  which  is  quite  a  lot  of 
scans,  this  implementation  only  requires  10  bits  of  local  site  data  (6  bits  for 
the  momentum  states  and  4  bits  for  messenger  states)  which  means  that  only 
1  kilobyte  is  needed  to  store  a  long-range  rule.  Since  the  CAM-8  uses  a  16-bit 
word,  there  still  remains  6  bits  of  unused  local  data.  We  use  these  remaining 
bits  to  hold  a  table  look-up  address.  That  is,  since  the  size  of  the  CAM-8 
look-up  table  static  random  access  memory  (SRAM)  is  64k  bytes,  and  our 
long-range  rule  only  requires  Ik  bytes,  we  can  store  up  to  64  long-range  rules 
into  CAM-8’s  SRAM  memory.  Since  our  molecular  dynamics  algorithm  de¬ 
composes  into  61  applications  of  the  long-range  rules,  it  is  now  clear  why 
we  have  chosen  to  use  up  to  eight  ranges.  Although  the  description  of  our 
implementation  may  sound  complicated,  in  fact  from  a  software  development 
point-of-view  it  was  the  most  direct  and  most  simple.  We  have  traded  off 
time  to  save  memory.  Yet  the  well  known  principle  of  computer  science  that 
one  can  save  much  time  at  the  expense  of  using  more  memory-  applies  to 
our  algorithm.  So  optimizations  of  our  algorithm  can  be  made,  particularly 
concerning  trading  off  an  increase  of  local  site  data  for  a  decrease  in  the  num¬ 
ber  of  needed  scans.  Clearly,  the  molecular  dynamics  algorithms  would  be 
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significantly  sped  up  if  thet  were  implemented  say  on  a  64-bit  architecture. 
Of  course  in  this  case,  computation  by  table  look-up  would  be  inappropri¬ 
ate.  However  by  making  use  of  lattice  isometries  and  rule  conjugation,  the 
necessary  logic  is  actually  quite  small,  as  evidenced  for  example  by  Eq.  (37) 
or  (52).  Therefore,  a  programmable  logic  method  of  computation  should  be 
better  than  the  table  look-up  method  of  computation  currently  in  use  in  the 
CAM-8  for  this  kind  of  lattice-gas  algorithm. 
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A  Computer  Implementation 


To  clearly  illustrate  the  method  of  long-range  lattice-gases,  we  have  provided 
a  computer  implementation  of  this  method  for  the  prototype  lattice-gas  ma¬ 
chine,  the  CAM-8.  Further  information  concerning  the  CAM-8  prototype  has 
been  provided  by  Margolus  [3].  The  language  we  use  to  implement  the  model 
is  called  “CAMForth”  and  was  developed  at  the  MIT  Laboratory  for  Com¬ 
puter  Science  as  a  superset  of  the  Forth  language  provided  on  SUN  SPARC- 
stations.  Local  computation  on  the  CAM-8  parallel  computer  is  achieved 
using  look-up  tables  stored  in  fast  static  ram  chips.  Since  the  CAM-8  data 
word  is  16  bits  wide,  a  lookup  table  contains  64A:  entries.  These  tables  are 
usually  created  off-line  and  stored  on  disk  as  binary  encoded  files,  where  the 
address  of  an  individual  entr}''  is  the  input  word,  and  the  value  of  that  entry 
is  the  output  word.  We  have  provided  some  C  language  subroutines  that  we 
use  to  create  the  look-up  tables  for  our  implementation.  On  the  eight-module 
CAM-8  prototype,  a  512  x  512  two-dimensional  simulation  runs  at  almost 
400  frames  per  second.  The  interface  dynamics,  that  is,  coalescence  of  drops, 
for  the  type  of  simple  liquid-gas  model  illustrated  here,  are  clearly  observable 
in  real-time  on  the  CAM-8’s  frame-buffer  display  or  within  an  X-window  on 
the  host  workstation. 
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A.l  Implementing  a  Lattice- Gas  in  the  CAMForth 
Language 

The  following  is  a  simple  example  of  a  long-range  lattice-gas  implemented 
in  the  CAMForth  language.  This  is  an  Appert-type  minimal  lattice-gas  [11] 
with  an  interaction  range  of  8.  It  is  implemented  here  with  periodic  bound¬ 
ary  conditions  in  two  dimensions  on  the  triangular  lattice.  The  example  code 
given  below  uses  a  512  by  512  lattice  size.  We  use  a  rhomboidal  mapping 
of  the  triangular  lattice  onto  the  square  lattice  by  using  lattice  directions: 
±f,  ±y,  and  ±{x  -1-  y).  Therefore,  on  the  CAM-8’s  frame-buffer  displaj^ 
the  representation  of  the  lattice-gas  on  screen  is  sheared.  Although  it  is 
possible  to  avoid  this,  for  simplicity  of  coding,  I  have  presented  this  exam¬ 
ple  which  is  the  most  direct  memor}^  mapping.  Our  implementation  uses  a 
single  messenger  bit,  or  "photon”,  for  each  particle  momentum  state.  There¬ 
fore,  our  implementation  uses  a  total  of  12  bits  of  local  site  data:  6  particle 
momentum  states  plus  6  photon  states.  Since  we  use  six  photon  states  the 
implementation  actually  does  simultaneous  updating  of  all  lattice  partitions. 

All  the  dynamics  are  encoded  into  two  16-bit  look-up  tables.  After  the 
particle  bits  are  kicked,  the  first  look-up  table,  Irlg.emit.lut,  actually  performs 
two  functions:  (1)  it  does  the  local  collisions;  and  (2)  it  then  emits  photons 
that  are  to  be  used  in  the  long-range  interaction  step.  After  the  photon  bits 
are  kicked,  the  second  look-up  table,  Mg. absorb,  lut,  performs  local  momen¬ 
tum  changes  by  using  the  kicked  photons,  and  in  this  way  locally  computes 
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the  long-range  interaction.  Any  interaction  range  can  be  chosen,  and  in  our 
code  given  below  is  set  in  the  constant  r.  If  a  negative  interaction  range  is 
set,  the  model  still  runs,  however  the  interaction  will  be  repulsive  instead 
of  being  attractive.  Therefore,  a  kind  of  parity  operation  with  respect  to 
the  interaction  range  causes  a  conjugation  of  the  interaction  operation  that 
results  in  flipping  the  polarity  of  the  interparticle  force. 
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\  *************  INITIALIZATION,  SITE  DATA,  AND  DATA  MOVEMENT 
new-experiment  512  by  512  space 

0  0  ==  pO  1  1  ==  pi  2  2  ==  p2  3  3  ==  p3  4  4  ==  p4  5  5  ==  p5 

6  6  ==  zO  7  7  ==  8  8  ==  z2  9  9  ==  z3  10  10  ==  z4  11  11  ==  z5 

:  p“kicks  kick  pO  field  -lx  0  y 
pi  field  0  X  -1  y 
p2  field  1  X  -1  y 
p3  field  lx  0  y 
p4  field  Ox  1  y 
p5  field  “lx  1  y 

8  constant  r 
r  negate  constant  -r 

:  z-kicks  kick  zO  field  -r  x  0  y 
zl  field  0  X  -r  y 
z2  field  r  X  -r  y 
z3  field  r  X  0  y 
z4  field  Ox  r  y 
z5  field  -r  X  r  y 

\  **********4:**  PREPARE  LOOK~UP  TABLES  AND  DEFINE  STEP 

create-lut  collision. tab  ""  Irlg . emit . lut  collision. tab  load-buffer 

create-lut  interaction. tab  ""  Irlg . absorb . lut  interact ion. tab  load-buffer 

:  prepare-tables 

lut-data  collision. tab  switch-luts 
lut-data  interaction. tab  switch-luts 
step 

;  this  is  when-starting 

define-step  stepx 

site-src  lut  lut-src  site 

p-kicks  run  new-table  z-kicks  run  new-table 
end-step  this  is  update-step 

\  *************  COLORMAP,  DISPLAY  TABLE,  AND  INITIAL  PATTERN 
""  Irlg. pal  palette  load-buffer  palette>display 
""  Irlg.dtab  display-table  load-buffer 
""  Irlg. pat  file>cam  xvds 
show-function  xvds 
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This  CAMForth  example  code,  Iglr.fth,  loads  several  additional  binary  en¬ 
coded  files  (to  set  up  for  frame-buffer  visualization,  data  interaction  through 
look-up  table,  and  initial  conditions)  to  run  on  the  CAM-8  prototype  or  the 
CAM-8  simulator.  These  additional  files  are  the  following: 

•  Irlg.pal  -  Color  palette  used  for  display 

•  Irlg.dtab  -  Look-up  table  to  render  the  density  field 

•  Irlg.pat  -  Initial  random  pattern  at  30 

•  Irlg.  emit,  lut  -  Local  collisions  and  photon  emission 

•  Irlg.absorb.lut  -  Long-range  interaction  table 


A.2  Creating  the  Look-up  Tables  in  the  C  Language 

I  have  written  a  short  subroutine  in  C  language  code  to  generate  the  look¬ 
up  tables  for  this  simple  example.  will  use  the  following  notation  for 
the  number  variables.  The  particle  momentum  states  are  n{a)  =  Ua  and 
_n(a)  =  l-na  and  np{a)  =  <  and  .np{a)  =  1  -  <•  Similarly  for  the  photon 
states,  z{a)  =  and  _z(o)  =  I  -  and  zp{a)  =  <  and  .zp[a)  =  1-4- 
The  number  variables  for  the  particles  and  photons  are  then  coded  as  global 
data  in  the  C  language  as  follows: 

#define  BIT(x,y)  ((y»x)%2) 

#defiiie  CAM^WORD  unsigned  short 
#define  B  6 

char  bit  [B]  ,  bitp  [B]  ; 
char  zbit  [B]  ,  zbitp  [B]  ; 

#define  n(a)  bit  [(a+B)7oB] 
tdefine  _n(a)  (l-n(a)) 

#define  np(a)  bitp  [(a+B)7oB] 

#define  _np(a)  (l-np(a)) 

#define  z(a)  zbit  [(a+B)7oB] 

#define  _z(a)  (l-z(a)) 


#define  zp(a)  zbitp  [(a+B)7oB] 
#define  .zp(a)  (l-zpCa)) 
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It  is  possible  to  express  the  long-range  lattice-gas  dynamics  as  a  logical 
function.  Furthermore,  it  is  possible  to  write  a  C  language  subroutine  as  a 
direct  implementation  of  this  logical  function  that  accepts  a  CAM  input  word 
and  computes  the  correct  CAM  output  word  returned  by  the  subroutine. 
The  first  look-up  table,  Irlg.emit.lut,  is  generated  by  calling  the  following  C 
language  routine  for  all  possible  input  configurations; 

CAM.WORD  lrlg_emit_lut  (CAM.WORD  input) 

{ 

char  a  ,  C  ; 

CAM_W0RD  output=0  ; 

/*  FHP  local  collisions  with  z-emission  */ 

/*  Convert  input  CAM  word  to  input  boolean  variables  */ 
for  (  a=0  ;  a<B  ;  a++)  n(a)  =  BIT(a,  input)  ; 

for  (  a=0  ;  a<B  ;  a++) 

C  =  n(a+2)  *  n(a+5)  *  _n(a+0)  *  _n(a+l)  *  _n(a+3)  *  _n(a+4) 

-  n(a+0)  *  n(a+3)  *  _n(a+l)  *  _n(a+2)  *  _n(a+4)  *  _n(a+5) 

+  n(a+l)  *  n(a+3)  *  n(a+5)  *  _n(a+0)  *  _n(a+2)  *  _n(a+4) 

-  _n(a+l)  *  _n(a+3)  *  _n(a+5)  *  n(a+0)  *  n(a+2)  *  n(a+4)  ; 

np(a)  =  n(a)  +  C  ; 

} 

/*  z-emission  */ 

for  (  a=0  ;  a<B  ;  a++)  z(a)  =  np(a+B/2)  *  _np(a)  ; 

/*  Convert  output  boolean  variables  to  output  CAM  word  */ 
for  (  a=0  ;  a<  B  ;  a++)  output  +=  np(a)«a  ; 
for  (  a=B  ;  a<2+B  ;  a++)  output  +=  z(a)«a  ; 


return  output  ; 

} 

The  low  six  bits  hold  the  particle  states,  bits  0  to  5,  and  the  high  six  bits  hold 
the  photon  states,  bits  6  to  11.  The  local  collision  operator  contains  the  FHP 
2-body  and  3-body  collisions  [1].  The  second  look-up  table,  Irlg.eTnit.lut,  is 
generated  by  calling  the  following  C  language  routine  for  all  possible  input 
configurations 
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CAM^WORD  Irlg^absorb.lut  (CAM^WORD  input) 

int  a  ; 

CAM^WDRD  oiitput=0  ; 
char  C  ; 


/*  Convert  input  CAM  word  to  input  boolean  variables  */ 
for  (  a=0  ;  a<B  ;  a++)  np(a)  =  BIT(a,  input)  ; 

/*  Condition  for  outgoing  photons  */ 

for  (  a=0  ;  a<B  ;  a++)  z(a)  =  np(a+B/2)  *  ^np(a)  ; 

/*  Streamed  incoming  photons  */ 

for  (  a=B  ;  a<2*B  ;  a++)  zpCa)  =  BIT(a,  input)  ; 

/*  Nonlocal  interactions,  photon  absorption  */ 
for  (  a=0  ;  a<B  ;  a++) 

C  =  zp(a+B/2)  *  z(a)  -  zp(a)  *  z(a+B/2)  ; 

n(a)  =  np(a)  +  C  ; 

} 

/*  Convert  output  boolean  variables  to  output  CAM  word  */ 
for  (  a=0  ;  a<B  ;  a++)  output  +=  n(a)<<a  ; 

return  output  ; 

} 

We  have  provided  the  algorithm  detail  necessary  for  one  to  understand  the 
essential  computational  structure  of  a  long-range  lattice-gas  model.  The 
example  presented  here  of  a  liquid-gas  system  was  coded  in  the  CAMForth 
language  and  the  C  language.  Code  to  generate  display  tables  and  color 
palettes  was  not  presented  here  since  that  is  not  essential  for  understanding 
our  lattice-gas  implementation  and  is  CAM-8  specific. 
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